La práctica se centra en el Parque Natural de Punta Entinas del Sabinar, una pequeña reseña del lugar la podemos encontrar en el siguiente enlace:
El proposito de esta práctica es hacer una fusion con las imagenes de Landsat. Por un lado, las imagenes mutiespectrales y por el otro la pancromática de resolución menor (15 metros).
library(raster)
## Loading required package: sp
library(terra)
## terra 1.7.3
library(rgdal)
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-4, (SVN revision 1196)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/javiv/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:\Program Files\PostgreSQL\13\share\contrib\postgis-3.0\proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
##
## Attaching package: 'rgdal'
## The following object is masked from 'package:terra':
##
## project
setwd("C:/MASTER_GOFOREST/2023/SENSORES/tarea_isable_2")
LST <- list.files(".", recursive = TRUE, full.names= TRUE, pattern = "B[12345678].TIF")
LST
## [1] "./LC08_L1TP_200034_20230228_20230228_02_RT_B1.TIF"
## [2] "./LC08_L1TP_200034_20230228_20230228_02_RT_B2.TIF"
## [3] "./LC08_L1TP_200034_20230228_20230228_02_RT_B3.TIF"
## [4] "./LC08_L1TP_200034_20230228_20230228_02_RT_B4.TIF"
## [5] "./LC08_L1TP_200034_20230228_20230228_02_RT_B5.TIF"
## [6] "./LC08_L1TP_200034_20230228_20230228_02_RT_B6.TIF"
## [7] "./LC08_L1TP_200034_20230228_20230228_02_RT_B7.TIF"
## [8] "./LC08_L1TP_200034_20230228_20230228_02_RT_B8.TIF"
LST <- lapply(1:length(LST), function (x) {raster(LST[x])})
#cargamos en un stack las bandas multiespectrales y por otro lado
#la banda Pancromatica
B8<-LST[[8]]
LST<-LST[1:7]
LST_stack <- stack(LST)
# we do a clip centered in our region of interest in this case
#is the Natural Park "Punta entinas del Sabinar" located in El Ejido, Almería.
PES<-readOGR("./PES.shp")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile
## Source: "C:\MASTER_GOFOREST\2023\SENSORES\tarea_isable_2\PES.shp", layer: "PES"
## with 1 features
## It has 6 fields
En este caso se trata del Paque Natural de Punta Entinas del Sabinar, unicado entre los municipios de El Ejido y Roquetas de Mar.
#we could not do a stack of all of our images because Pancromatic band
#has less scale or size of pixel compared with the others, the multiespectral images.
PES_stack_crop <- terra::crop(LST_stack, PES)
plot(PES_stack_crop)
B8<-terra::crop(B8,PES)
plot(B8)
# 1- DATA REESCALE
PES_15<-resample(PES_stack_crop,B8)
dim(PES_15)
## [1] 375 851 7
dim(B8)
## [1] 375 851 1
Procedemos a calcular las estadisticas media y desviacion tipica de cada una de las bandas multiespectral.
#Calculate mean and standar desviation of each band of multiespectral images
#extract values of each band of image
sts <- data.frame(B_mean = numeric(),B_sd = numeric())
for (i in 1:ncol(PES_15@data@values)){
sts[i,1] <-mean(PES_15@data@values[,i])
sts[i,2] <-sd(PES_15@data@values[,i])
rownames(sts)[i]<-paste0("B",i, sep="")
}
#añadimos las estadisticas mean y sd de la banda pancromatica a tabla de estadisticas
sts[c("B8"),1]<-c(mean(B8[]))
sts[c("B8"),2]<-c(sd(B8[]))
#Transformamos la informacion de los píxeles de la imagen numericamente con la Fórmula de Torus.
# Calculate the parameters of the transformation image: Torus's Method
#we do the same loop to calculate ai and bi from each spectral band
# ai = sd_Bi/sd_B8_PAN
# bi = mean_Bi-((sd_Bi/sd_B8_PAN)*mean_B8_PAN)
# And finally:
#PAN_Bi = ai * PAN +bi
# WE CALCULATE
#define a data.frame with a and b like empty comlumns
rm(param)
## Warning in rm(param): objeto 'param' no encontrado
param<-data.frame(a= numeric(),
b= numeric())
#calculate each index for each spectral band
for (i in 1:ncol(PES_15@data@values)){
param[i,1] <-sts$B_mean[i]/sts$B_mean[8]
param[i,2] <-sts$B_mean[i]-((sts$B_sd[i]/sts$B_sd[8])*sts$B_mean[8])
rownames(param)[i]<-paste0("B",i, sep="")
}
#PAN_B1 = a1 * B8 +b1
rm(PAN_stk)
## Warning in rm(PAN_stk): objeto 'PAN_stk' no encontrado
PAN_stk<-list()
for (i in 1:ncol(PES_15@data@values)){
PAN_stk[i] <-param$a[i]*B8 +param$b[i]
names(PAN_stk[[i]])<-paste0("B",i, sep="")
}
## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated
PAN_stk<-stack(PAN_stk)
plot(PAN_stk)
#Define the filter (kernel1)
# Apply the filter that will smooth each image
kernel1 <- matrix(c(0.0039060000, 0.015625000, 0.023438000, 0.015625000, 0.0039060000,
0.015625000, 0.062500000, 0.093750000, 0.062500000, 0.015625000,
0.023438000, 0.093750000, 0.14062500, 0.093750000, 0.023438000,
0.015625000, 0.062500000, 0.093750000, 0.062500000, 0.015625000,
0.0039060000, 0.015625000, 0.023438000, 0.015625000, 0.0039060000),
nrow = 5)
print(kernel1)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.003906 0.015625 0.023438 0.015625 0.003906
## [2,] 0.015625 0.062500 0.093750 0.062500 0.015625
## [3,] 0.023438 0.093750 0.140625 0.093750 0.023438
## [4,] 0.015625 0.062500 0.093750 0.062500 0.015625
## [5,] 0.003906 0.015625 0.023438 0.015625 0.003906
#Apply the filter to all multiespectral images transformed
filtered<-list()
for (i in 1:ncol(PES_15@data@values)){
filtered[i] <-focal(PAN_stk[[i]],kernel1)
names(filtered[[i]])<-paste0("B",i, sep="")
}
## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated
filtered_stack<-stack(filtered)
plot(filtered_stack)
#When we alpply the low filter o the images next step is subtract this images with originals images
Detail<-list()
for (i in 1:ncol(PES_15@data@values)){
Detail[i] <-PAN_stk[[i]]-filtered[[i]]
names(Detail[[i]])<-paste0("DetB",i,sep="")
}
## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated
Detail_stk<-stack(Detail)
plot(Detail_stk)
### SUMA DE BANDAS: DETALLE + ORIGINAL
#once we applied the "filtro por lo alto" we sum the bands
fusion<-list()
for (i in 1:ncol(PES_15@data@values)){
fusion[i] <-PES_15[[i]]+Detail_stk[[i]]
names(fusion[[i]])<-paste0("FusB",i,sep="")
}
## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated
## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated
fusion_stk<-stack(fusion)
plot(fusion_stk)
Representamos las bandas del infrrarrojo, del rojo y el verde.
plotRGB(fusion_stk, r=5, g=4, b=3, scale=maxValue(fusion_stk[[5]]), stretch= "hist")
#writeRaster(fusion_stk,"./stack_fus_automatic.tiff")
#una vez que esta echa la funsion de las imagenes puedo hacer una clasificaion no supervisada y
#comparar si existen variaciones en la clasificacion gracias a la fusion de la imagen
b<-fusion_stk
b[is.na(b)]<-0
kMeansResult <- kmeans(scale(b[]), centers=6)
kMeansResult$centers
## FusB1 FusB2 FusB3 FusB4 FusB5 FusB6 FusB7
## 1 0.3163838 0.3307411 0.3548348 0.37616767 0.4628357 0.2950999 0.1384995
## 2 -0.6618898 -0.6953083 -0.6977723 -0.62983171 -0.3802437 -0.2323143 -0.1768918
## 3 1.7718918 1.7941766 1.7616567 1.71606633 1.4940061 1.4701645 1.4996942
## 4 1.0441380 1.0686035 1.0752250 1.07383480 1.0003096 0.8885167 0.8020898
## 5 -1.0160851 -1.0339415 -1.0872439 -1.18219236 -1.3731142 -1.3986488 -1.3285032
## 6 -0.3302744 -0.3275334 -0.2247280 -0.08284935 0.1681415 0.5966643 0.7992552
#creamos el raster vacio para luego rellenarlo
result <- raster(b[[1]])
#rellenamos el raster creado con los datos del cluster del k-means
result <- setValues(result, kMeansResult$cluster)
#writeRaster(result,"./stack_fus_kmeans_1.tiff")
plot(result, col=c("red","black","yellow","green","blue","white"))
#ahora probamos la clasificacion con las imagenes normales sin modificar
b<-PES_stack_crop
b[is.na(b)]<-0
kMeansResult <- kmeans(scale(b[]), centers=6)
kMeansResult$centers
## LC08_L1TP_200034_20230228_20230228_02_RT_B1
## 1 0.3236445
## 2 1.8669920
## 3 -0.8137024
## 4 1.0911504
## 5 -1.0492191
## 6 -0.4410001
## LC08_L1TP_200034_20230228_20230228_02_RT_B2
## 1 0.3272035
## 2 1.8567103
## 3 -0.8242341
## 4 1.0920240
## 5 -1.0532609
## 6 -0.4167034
## LC08_L1TP_200034_20230228_20230228_02_RT_B3
## 1 0.3457095
## 2 1.7922517
## 3 -0.8069011
## 4 1.0807066
## 5 -1.1028491
## 6 -0.2770841
## LC08_L1TP_200034_20230228_20230228_02_RT_B4
## 1 0.3614099
## 2 1.7198005
## 3 -0.7065811
## 4 1.0632366
## 5 -1.1984324
## 6 -0.1051973
## LC08_L1TP_200034_20230228_20230228_02_RT_B5
## 1 0.4536012
## 2 1.4626346
## 3 -0.4065128
## 4 0.9707657
## 5 -1.3977319
## 6 0.1792124
## LC08_L1TP_200034_20230228_20230228_02_RT_B6
## 1 0.2807386
## 2 1.4214960
## 3 -0.2322147
## 4 0.8637451
## 5 -1.4336181
## 6 0.6577859
## LC08_L1TP_200034_20230228_20230228_02_RT_B7
## 1 0.1151378
## 2 1.4539015
## 3 -0.1719946
## 4 0.7788549
## 5 -1.3668828
## 6 0.9017684
#creamos el raster vacio para luego rellenarlo
result <- raster(b[[1]])
#rellenamos el raster creado con los datos del cluster del k-means
result <- setValues(result, kMeansResult$cluster)
#writeRaster(result,"./landsat_kmeans.tiff")
plot(result, col=c("red","black","yellow","green","blue","white"))